Magnetic resonance imaging method

ABSTRACT

In SENSitivity Encoding (SENSE), the reconstructed images are susceptible to amplified noise and/or artifacts if the underlying matrix inversion procedure is ill-conditioned. In this work, we propose to firstly apply the conventional SENSE algorithm to obtain an initial estimate. This initial estimate undergoes filtering to improve the signal-to-noise ratio. Then, it is fed back to the reconstruction as a reference image to estimate the amount of aliasing that may arise from regularization. We derive the optimal regularized solution that minimizes the weighted sum of artifact and noise power.

The present invention relates to a magnetic resonance imaging methodwherein undersampled magnetic resonance signals are acquired by areceiver antennae system having a spatial sensitivity profile and theimage being reconstructed from the undersampled magnetic resonancesignals and the spatial sensitivity profile.

In the method of undersampled acquisition known as SENSitivity Encoding(SENSE), as described by Pruessmann K P et al. in MRM 1999; 42:952-962,the reconstructed images are susceptible to amplified noise and/orartifacts if the underlying matrix inversion procedure isill-conditioned.

It is therefore an object of the present invention to improve the abovementioned method of undersampled acquisition by reducing the amount ofnoise in the final image in order to handle scenarios with a lowsignal-to-noise ratio (SNR) and/or a high reduction factor R.

This object is achieved by the magnetic resonance imaging methodaccording to claim 1, whereas the reconstruction of the image isprovided by a first step, in which image is reconstructed on the basisof reconstruction matrices according to a parallel imaging like SENSE,thereinafter the so reconstructed image is subject to a filteringoperation, which provides a post-processed image, which is used to alterthe reconstruction matrices, and by a second step, in which the finalimage is reconstructed on the basis of the altered reconstructionmatrices. Using this scheme, the reconstruction from the second step isoptimized with respect to minimizing noise and aliasing artifacts. Thefurther objects are achieved by the magnetic resonance imaging systemaccording to claim 8 and the computer program product according to claim9.

The method according to the present invention has the advantage, thatthe amount of noise artifacts in the image can be reduced without anyinfluence on the sampling rate, i.e. the reduction factor R.

These and other aspects of the invention will be elaborated withreference to the preferred implementations as defined in the dependentclaims. In the following description an exemplified embodiment of theinvention is described with respect to the accompanying drawings. Itshows

FIG. 1 a diagram of the acceleration factor R versus the normalized RMSerror (left) and an reconstructed image with SENSE only and withfeedback regularization (right), and

FIG. 2 diagrammatically a magnetic resonance imaging system in which theinvention is used.

In the present description a multiple of receiver antenna or coils areused. However, it is also possible to implement the SENSE method with asingle receiving coil or antenna at different receiving positions.

Basic Principles

At the heart of the SENSitivity Encoding (SENSE), a parallel imaging(PI) method, as described in Pruessmann K P, et al. Magn Reson Med42:952-962, 1999, lies a series of matrix inversions that determine theunaliased image voxels v from the measured k-space data a This linearsystem can be represented as:v=(L ⁻¹ S)⁺ L ⁻¹ a=F a   [1]where S denotes the so-called sensitivity matrix (1). L is the “squareroot” (e.g. by Cholesky decomposition) of the noise correlation matrix Ψ(2) (i.e. Ψ=L L^(H)). Superscript+denotes a (regularized)pseudo-inverse. Eq. [1] is equivalent to the original SENSE formulationas described in Pruessmann K P, et al. Magn Reson Med 42:952-962, 1999and in Pruessmann K P, et al. Magn Reson Med 46:638-651, 2001.

If the matrix product (L⁻¹ S) is ill-conditioned, v is sensitive toperturbations (cf. Golub G H, Van Loan C F. Matrix computations. 3^(rd)ed. Johns Hopkins University Press, 1996) on the right hand side of Eq.[1], including measurement noise and inaccuracy of the sensitivity maps.A wide variety of regularization approaches exist to improve theconditioning of the inversion procedure (see e.g. Hansen P C. NumericalAlgorithms 6:1-35, 1994). In all cases, accuracy of the inversion istraded off to gain stability. In the present work, we propose a two-passprocedure to estimate this trade-off quantitatively.

PRACTICAL EXAMPLE OF THE INVENTION

In the first pass of the proposed method, the conventional SENSEalgorithm is applied using only truncated singular value decomposition(SVD) to avoid obvious noise amplification (cutoff at conditionnumber >100). This generates an initial estimate {circumflex over (v)},which undergoes median filtering to improve the signal-to-noise ratio.In the second pass, the regularized reconstruction matrix F isdetermined as the solution that minimizes the following weighted sum:α(Noise Power)+(Artifact Power)=α∥F L∥ _(Frob) ²+∥(FS−I)diag({circumflex over (v)})∥_(Frob) ²   [2]where α denote the weight given to the noise term relative to theartifact term (arbitrarily set to 1 in this work); ∥•∥_(Frob) denotesthe Frobenius norm. The first term of Eq. [1] estimates the noise powerof the reconstructed voxels, while the second term estimates theartifact power resulting from regularization assuming that the truevoxel intensities are given by {circumflex over (v)}. Regardless of theregularization strategy used (e.g. diagonal loading, truncated SVD,damped SVD, etc.) the optimal reconstruction matrix F_(opt) thatminimizes Eq. [2] can be determined analytically, and it has severalmathematically equivalent forms, including: $\begin{matrix}\begin{matrix}{F_{opt} = {\left( {{S^{H}\psi^{- 1}S} + {\alpha\quad{{diag}\left( {\hat{v}}^{2} \right)}^{- 1}}} \right)^{- 1}S^{H}\psi^{- 1}}} \\{= {{{diag}\left( {\hat{v}}^{2} \right)}{S^{H}\left( {{S\quad{{diag}\left( {\hat{v}}^{2} \right)}S^{H}} + {\alpha\quad\psi}} \right)}^{- 1}}}\end{matrix} & \left\lbrack {{3\quad a},b} \right\rbrack\end{matrix}$

For α=1, these expressions become equivalent to those previously derived(8-9). An interesting observation is that the matrix product Sdiag({circumflex over (v)}) is equal to the sensitivity maps multipliedby image estimate {circumflex over (v)}. This has been referred to asthe “in vivo sensitivities” (cf. Wang J et al. Workshop on Parallel MRImaging Basics and Clinical Applications. 89, 2001, and Sodickson D K.Magn Reson Med. 44:243-251, 2000). Thus, Eq. [3b] can be rewritten asfollows, with S_(in vivo)=S diag({circumflex over (v)}):F _(opt)=diag({circumflex over (v)})S _(in vivo) ^(H)(S _(in vivo) S_(in vivo) ^(H)+α↓)⁻¹   [4]

In principle, the use of in vivo sensitivities has no effect on thereconstruction. In practice however, the in vivo sensitivities aretypically acquired using the center of k-space (compare McKenzie C A, etal. Workshop on Parallel MR Imaging Basics and Clinical Applications.88, 2001) or a separate low-resolution reference. Thus, the in vivosensitivities are convolved with a low-pass point spread function. Thisapproximation can be regarded as a modeling error in Eq. [1]. Themaximum error amplification is bounded by the condition number ofF_(opt) (see Golub G H, Van Loan C F. Matrix computations. 3 ed.Baltimore: Johns Hopkins University Press, 1996.); while the minimumerror is bounded by the reconstruction error from actually usingaccurate high-resolution in vivo sensitivities as S_(in vivo).

Simulations were performed using a cardiac image (see e.g. Weiger M, etal. Magn Reson Med 43:177-184, 2000), with a six-element coil arrayplaced around the body (see e.g. Weiger M, et al. Magn Reson Med45:495-504, 2001), and a signal-to-noise ratio of 10. Root-mean-square(RMS) reconstruction error was determined as a function of theacceleration factor (R) along the phase-encoding direction (left-right).

FIG. 1 shows that the RMS error improves at all acceleration factorswith regularization, including a marginal improvement even at R=1. Thisimprovement at R=1 is due to the feedback mechanism serving as a aself-consistency check. In general, the amount of improvement stronglydepends on the image contents, with larger improvements possible if thealiased voxels exhibit high contrasts. On the other hand, if the entirefield-of-view has approximately uniform intensity, the improvement isnegligible, as would be expected. The inset images show thereconstruction results with and without feedback regularization atR=4.5.

CONCLUSION

In the present invention, we present a feedback framework forregularized reconstruction. We exploit the fact that neighboring voxelsare highly correlated. As a result, filtering applied to the first-passreconstruction can be used to obtain a high signal-to-noise imageestimate, which can be used to estimate the potential amount ofartifacts. In the case of dynamic imaging, temporal correlations (asdescribed in Wang J et al. Workshop on Parallel MR Imaging Basics andClinical Applications. 89, 2001) or joint spatiotemporal correlationsmay also be used to obtain the image estimate. For a given estimate, weapplied median filtering to improve the image quality, but a widevariety of other filtering methods can be used as well, includinganisotropic diffusion (see Gerig G, et al. IEEE Trans Med Imaging11:221-232, 1992) and statistical approaches. Finally, thenoise-versus-artifact tradeoff can also be evaluated in a number ofmanners (see e.g. Hansen P C. Numerical Algorithms 6:1-35, 1994).However, for any regularized reconstruction that minimizes theexpression in Eq. [2], the reconstruction formula in Eqs. [3a, 3b, 4]represent the optimal, regardless of the regularization strategy used.

In the above presented method a number of filtering methods can be used,such as low-pass filtering, median filtering, statistical filtering,anisotropic filtering and wavelet filtering. Low-pass filtering involvesblurring each voxel with its neighbours. Median filtering involvesreplacing the intensity of each voxel wiht the median of the voxelintensities within a neighbourhood. Statistical filtering involvescomparing the statistical properties of each voxel to those of noise,and discarding or attenuating those voxels that are similar to noise.Anisotropic filtering involves blurring each voxel with its neighbourswith the degree of blurring dependent on the degree of similaritybetween them. Wavelet filtering involves transforming an image fromgeometric space to wavelet space, which is spanned by a family ofwavelet functions. The filtering is then applied in wavelet space usingany of the above filtering methods. The filtered data areinverse-transformed back to geometric space.

FIG. 3 shows diagrammatically a magnetic resonance imaging System inwhich the invention is used.

The magnetic resonance imaging system includes a set of main coils 10whereby a steady, uniform magnetic field is generated. The main coilsare constructed, for example in such a manner that they enclose atunnel-shaped examination space. The patient to be examined is slid on atable into this tunnel-shaped examination space. The magnetic resonanceimaging system also includes a number of gradient coils 11, 12 wherebymagnetic fields exhibiting spatial variations, notably in the form oftemporary gradients in individual directions, are generated so as to besuperposed on the uniform magnetic field. The gradient coils 11, 12 areconnected to a controllable power supply unit 21. The gradient coils 11,12 are energized by application of an electric current by means of thepower supply unit 21. The strength, direction and duration of thegradients are controlled by control of the power supply unit. Themagnetic resonance imaging system also includes transmission andreceiving coils 13, 15 for generating RF excitation pulses and forpicking up the magnetic resonance signals, respectively. Thetransmission coil 13 is preferably constructed as a body coil whereby (apart of) the object to be examined can be enclosed. The body coil isusually arranged in the magnetic resonance imaging system in such amanner that the patient 30 to be examined, being arranged in themagnetic resonance imaging system, is enclosed by the body coil 13. Thebody coil 13 acts as a transmission aerial for the transmission of theRF excitation pulses and RF refocusing pulses. Preferably, the body coil13 involves a spatially uniform intensity distribution of thetransmitted RF pulses. The receiving coils 15 are preferably surfacecoils 15 which are arranged on or near the body of the patient 30 to beexamined. Such surface coils 15 have a high sensitivity for thereception of magnetic resonance signals which is also spatiallyinhomogeneous. This means that individual surface coils 15 are mainlysensitive for magnetic resonance signals originating from separatedirections, i.e. from separate parts in space of the body of the patientto be examined. The coil sensitivity profile represents the spatialsensitivity of the set of surface coils. The transmission coils, notablysurface coils, are connected to a demodulator 24 and the receivedmagnetic resonance signals (MS) are demodulated by means of thedemodulator 24. The demodulated magnetic resonance signals (DMS) areapplied to a reconstruction unit. The reconstruction unit reconstructsthe magnetic resonance image from the demodulated magnetic resonancesignals (D)MS) and on the basis of the coil sensitivity profile of theset of surface coils. The coil sensitivity profile has been measured inadvance and is stored, for example electronically, in a memory unitwhich is included in the reconstruction unit. The reconstruction unitderives one or more image signals from the demodulated magneticresonance signals (DMS), which image signals represent one or more,possibly successive magnetic resonance images. This means that thesignal levels of the image signal of such a magnetic resonance imagerepresent the brightness values of the relevant magnetic resonanceimage. The reconstruction unit 25 in practice is preferably constructedas a digital image processing unit 25 which is programmed so as toreconstruct the magnetic resonance image from the demodulated magneticresonance signals and on the basis of the coil sensitivity profile. Thedigital image processing unit 25 is notably programmed so as to executethe reconstruction in conformity with the so-called SENSE technique orthe so-called SMASH technique. The image signal from the reconstructionunit is applied to a monitor 26 so that the monitor can display theimage information of the magnetic resonance image (images). It is alsopossible to store the image signal in a buffer unit 27 while awaitingfurther processing, for example printing in the form of a hard copy.

In order to form a magnetic resonance image or a series of successivemagnetic resonance images of the patient to be examined, the body of thepatient is exposed to the magnetic field prevailing in the examinationspace. The steady, uniform magnetic field, i.e. the main field, orientsa small excess number of the spins in the body of the patient to beexamined in the direction of the main field. This generates a (small)net macroscopic magnetization in the body. These spins are, for examplenuclear spins such as of the hydrogen nuclei (protons), but electronspins may also be concerned. The magnetization is locally influenced byapplication of the gradient fields. For example, the gradient coils 12apply a selection gradient in order to select a more or less thin sliceof the body. Subsequently, the transmission coils apply the RFexcitation pulse to the examination space in which the part to be imagedof the patient to be examined is situated. The RF excitation pulseexcites the spins in the selected slice, i.e. the net magnetization thenperforms a precessional motion about the direction of the main field.During this operation those spins are excited which have a Larmorfrequency within the frequency band of the RF excitation pulse in themain field. However, it is also very well possible to excite the spinsin a part of the body which is much larger man such a thin slice; forexample, the spins can be excited in a three-dimensional part whichextends substantially in three directions in the body. After the RFexcitation, the spins slowly return to their initial state and themacroscopic magnetization returns to its (thermal) state of equilibrium.The relaxing spins then emit magnetic resonance signals. Because of theapplication of a read-out gradient and a phase encoding gradient, themagnetic resonance signals have a plurality of frequency componentswhich encode the spatial positions in, for example the selected slice.The k-space is scanned by the magnetic resonance signals by applicationof the read-out gradients and the phase encoding gradients. According tothe invention, the application of notably the phase encoding gradientsresults in the sub-sampling of the k-space, relative to a predeterminedspatial resolution of the magnetic resonance image. For example, anumber of lines which is too small for the predetermined resolution ofthe magnetic resonance image, for example only half the number of lines,is scanned in the k-space.

1. Magnetic resonance imaging method for forming, wherein undersampledmagnetic resonance signals are acquired by at least one receiver antennahaving a plurality of receiver antenna positions, each with a spatialsensitivity profile, the image being reconstructed from the undersampledmagnetic resonance signals and the spatial sensitivity profiles, whereasthe reconstruction of the image is provided by a first step, in whichthe image is reconstructed on the basis of reconstruction matricesaccording to a parallel imaging method, thereinafter the soreconstructed image is subject to a filtering operation, which providesa post-processed image, which is used to alter the reconstructionmatrices, and by a second step, in which the final image isreconstructed on the basis of the altered reconstruction matrices: 2.Magnetic resonance method according to claim 1, wherein the filtering isa median filtering method.
 3. Magnetic resonance method according toclaim 1, wherein the filtering is a wavelet filtering.
 4. Magneticresonance method according to claim 1, wherein the filtering is alow-pass filtering.
 5. Magnetic resonance method according to claim 1,wherein the filtering is performed by anisotropic diffusion.
 6. Magneticresonance method according to claim 1, wherein the filtering isperformed statistically.
 7. Magnetic resonance method according to claim1, wherein the alteration of the reconstruction matrix is performedbased on adjusting the trade-off between the noise level and theartifact level.
 8. A magnetic resonance imaging system comprising astatic main magnet having a main magnetic field, at least one receiverantenna having a plurality of receiver antenna positions, means forapplying a read and other gradients, means for measuring MR signalsalong a predetermined trajectory containing a plurality of lines ink-space a receiver antenna system for acquiring undersampled MR signals,each receiver antenna position having a spatial sensitivity profile,means for reconstruction of the image in a first step on the basis ofrecconstruction matrices according to a parallel imaging method, meansfor filtering the so reconstructed image, which provides apost-processed image, means for altering the reconstruction matrices byuse of the post-processed images, means for reconstruction of the finalimage on the basis of the altered reconstruction matrices.
 9. A computerreadable media comprising instructions for controlling a computer systemto perform a method of forming a magnetic resonance image the methodcomprising: applying a read and other gradients, measuring MR signalsalong a predetermined trajectory containing a plurality of lines ink-space acquiring undersampled MR signals from a receiver antennasystem, each receiver antenna position having a spatial sensitivityprofile, reconstruction of the image in a first step on the basis ofrecconstruction matrices according to a parallel imaging method,filtering the so reconstructed image, which provides a post-processedimage, altering the reconstruction matrices by use of the post-processedimages, reconstruction of the final image on the basis of the alteredreconstruction matrices.